def EGtest(DS_input):
"""
Tests if a given degree sequence DS_input can form a simple graph according to the Erdős–Gallai theorem.
Parameters:
- DS_input : the degree sequence of a graph.
Returns:
- 'success' if the sequence can form a simple graph, 'failure' otherwise.
"""
# Check if the sum of degrees is odd or if the sequence cannot form a simple graph
if np.sum(DS_input) % 2 == 1 or np.sum(DS_input) > len(DS_input) * (len(DS_input) - 1):
return 'failure'
# Remove zeros for proper delete operation, if necessary
if np.sum(DS_input) != 0 and 0 in DS_input:
DS_input = np.delete(DS_input, np.where(DS_input == 0))
# Sort in non-increasing order
DS_input = np.sort(DS_input)[::-1]
# Apply the Erdős–Gallai theorem
for index in range(len(DS_input)):
k = index + 1
tmp_DS_input = np.array([
k if DS_input[x] > k and x >= k else DS_input[x] for x in range(len(DS_input))
])
if sum(tmp_DS_input[0:k]) > k * (k - 1) + np.sum(tmp_DS_input[k:]):
return 'failure'
return 'success'
def rightmost_adj(DS_input):
"""
Identifies the rightmost adjacency set for a graph based on its degree sequence.
Parameters:
- DS_input : the degree sequence of a graph.
Returns:
- rightmost_adj_set : indices of the rightmost adjacency set.
"""
non_zero_degree_node_index = np.array([i for i, x in enumerate(DS_input) if x > 0])
if 0 in DS_input:
DS_input = np.delete(DS_input, np.where(DS_input == 0))
order_index = np.argsort(-DS_input) # Descending order index
DS_input = np.sort(DS_input)[::-1] # Sort in descending order
rightmost_adj_set = np.array([], dtype=int)
for non_leading_node in range(1, len(DS_input))[::-1]: # Iterate in reverse order
tmp_DS = np.copy(DS_input)
tmp_DS[0] -= 1 # Decrease the degree of the leading node
tmp_DS[non_leading_node] -= 1 # Decrease the degree of the non-leading node
# Prepare for Erdős–Gallai test
if tmp_DS[0] != 0:
DS_for_test = np.array([tmp_DS[i] - 1 if i <= tmp_DS[0] else tmp_DS[i] for i in range(1, len(tmp_DS))])
else:
DS_for_test = np.copy(tmp_DS)
if EGtest(DS_for_test) == 'success':
rightmost_adj_set = np.append(rightmost_adj_set, non_leading_node)
if tmp_DS[0] == 0:
break # Stop if the leading node's degree reaches 0
DS_input = np.copy(tmp_DS)
# Reorder indices to match original sequence
rightmost_adj_set = np.array(non_zero_degree_node_index[order_index[rightmost_adj_set]])
return rightmost_adj_set
def left_adj(current_adj_mat, DS_input_original, DS_input_current, rightmost_adj, leading_node_index):
"""
Determines the left adjacency set for the current graph.
Parameters:
- current_adj_mat : the current adjacency matrix.
- DS_input_original : the original degree sequence.
- DS_input_current : the current degree sequence after updates.
- rightmost_adj : indices of the rightmost adjacency set.
- leading_node_index : index of the leading node.
Returns:
- left_adj_set : indices of nodes forming the left adjacency set.
"""
leading_node_degree = DS_input_current[leading_node_index]
tmp_DS = np.copy(DS_input_current)
tmp_DS[leading_node_index] = 0
order_index = np.argsort(-tmp_DS)
tmp_DS = np.sort(tmp_DS)[::-1]
# remove node withh zero degree
non_zero_DS = np.delete(tmp_DS, np.where(tmp_DS == 0))
number_of_min = len(np.where(non_zero_DS == min(non_zero_DS))[0])
non_min_node_index = np.where(non_zero_DS != min(non_zero_DS))[0]
number_of_non_min = len(non_min_node_index)
# parameters for for-loop
tmp_left_adj_set = [[0]*leading_node_degree]
i_start = max(leading_node_degree - number_of_min, 0) # 0要檢查 因為python index 0 = R inedx 1
i_end = max(min(number_of_non_min, leading_node_degree), 0) + 1
# deal with duplicated structure
duplicate_marker = [mat.copy() for mat in current_adj_mat]
[duplicate_marker[i].append(DS_input_original[i]) for i in range(len(DS_input_original))]
duplicated_index = [0]*len(duplicate_marker)
for i in range(len(duplicate_marker)):
for j in range(len(duplicate_marker)):
if duplicate_marker[j] == duplicate_marker[i]:
duplicated_index[i] = j
duplicated_index = [duplicated_index[order_index[i]] for i in range(len(order_index))]
for i in range(i_start, i_end):
# first part (for non-min degree node)
if i == 1 and number_of_non_min == 1:
first_part = [list(non_min_node_index)]
elif i != 0:
first_part = [list(l) for l in combinations(non_min_node_index,i)]
duplicated_mat = [0]*len(first_part)
for j in range(len(first_part)):
duplicated_mat[j] = [duplicated_index[first_part[j][k]] for k in range(len(first_part[j]))]
unique_index = []
unique_value = []
for j in range(len(duplicated_mat))[::-1]:
x = duplicated_mat[j]
if x not in unique_value:
unique_value.append(x)
unique_index.append(j)
first_part = [first_part[m] for m in range(len(first_part)) if m in unique_index]
# second part (for min degree node)
if i != leading_node_degree:
min_degree_node = np.where(non_zero_DS == min(non_zero_DS))[0]
if len(min_degree_node) == 1:
second_part = [list(min_degree_node)]
else:
second_part = [list(l) for l in combinations(min_degree_node,leading_node_degree-i)]
duplicated_mat = [0]*len(second_part)
for j in range(len(second_part)):
duplicated_mat[j] = [duplicated_index[second_part[j][k]] for k in range(len(second_part[j]))]
unique_index = []
unique_value = []
for j in range(len(duplicated_mat))[::-1]:
x = duplicated_mat[j]
if x not in unique_value:
unique_value.append(x)
unique_index.append(j)
second_part = [second_part[m] for m in range(len(second_part)) if m in unique_index]
#combine first part & second part
if i == 0:
combine_two_part = second_part
elif i == leading_node_degree:
combine_two_part = first_part
else:
combine_two_part = [x + y for x in first_part for y in second_part]
tmp_left_adj_set = tmp_left_adj_set + combine_two_part
tmp_left_adj_set.remove(tmp_left_adj_set[0]) # delete offset
## check colex order
# calculate colex score (ex: colex_order_rightmost = [0 2 3] -> clex_score_rightmost = [2^0 2^2 2^3] = [1 4 8])
# ex: tmp_left_adj_set = [[1,2,3],[0,2,3]] -> score = 14, 13
# colex_order_rightmost = [0, 2, 3] -> score = 13
# then exclude [1,2,3]
mapping_index = [np.where(tmp_DS == tmp_DS[k])[0][0] for k in range(len(tmp_DS))]
colex_order_rightmost = [np.where(order_index == rightmost_adj[i])[0][0] for i in range(len(rightmost_adj))]
colex_score_rightmost = np.sum([2**(mapping_index[colex_order_rightmost[i]]) for i in range(len(colex_order_rightmost))])
colex_score_left = []
for i in range(len(tmp_left_adj_set)):
colex_score_i = np.sum([2**(mapping_index[tmp_left_adj_set[i][j]]) for j in range(len(colex_order_rightmost))])
colex_score_left.append(colex_score_i)
check_to_the_left = np.array(colex_score_left <= colex_score_rightmost, dtype=bool)
# filter
tmp_left_adj_set = np.asarray(tmp_left_adj_set)
left_adj_set = tmp_left_adj_set[check_to_the_left]
left_adj_set = [list(order_index[k]) for k in left_adj_set]
return left_adj_set
def connect_adj_set(leading_node_index, current_adj_mat, adj_set):
"""
Connects nodes in the adjacency set to the leading node, updating the adjacency matrix.
Parameters:
- leading_node_index : index of the leading node in the adjacency matrix.
- current_adj_mat : current adjacency matrix of the graph.
- adj_set : set of nodes to be connected to the leading node.
Returns:
- output_mat : updated adjacency matrix with new connections.
"""
output_mat = []
for ii in range(len(adj_set)):
tmp_mat = [row.copy() for row in current_adj_mat] # Create a deep copy of the current adjacency matrix
for jj in range(len(adj_set[0])):
# Update adjacency for both the leading node and the current node in the set
tmp_mat[leading_node_index][adj_set[ii][jj]] = tmp_mat[adj_set[ii][jj]][leading_node_index] = 1
output_mat.append([row.copy() for row in tmp_mat])
return output_mat
def net_gen(original_DS):
"""
Generates potential network structures from a given degree sequence.
Parameters:
- original_DS : the original degree sequence for network generation.
Returns:
- complete_adj_mat : generated potential network structures.
"""
sum_DS = np.sum(original_DS)
rows, cols = len(original_DS), len(original_DS)
incomplete_adj_mat = [[[0] * cols for _ in range(rows)]]
complete_adj_mat = []
max_complete_mat_count = 1000 # Limit for generated potential network structures of each degree sequence
while incomplete_adj_mat and len(complete_adj_mat) < max_complete_mat_count:
last_matrix = incomplete_adj_mat.pop() # Get the last matrix to work on
current_DS = original_DS - np.array([sum(row) for row in last_matrix]) # Update the degree sequence
if np.sum(current_DS != 0) > 1: # Ensure there's more than one node left to connect
leading_node = np.argmax(current_DS) # Find the node with the highest degree
rightmost_adj_set = rightmost_adj(current_DS) # Find rightmost adjacency set
# Generate left adjacency sets
left_adj_set = left_adj(
current_adj_mat=last_matrix,
DS_input_original=original_DS,
DS_input_current=current_DS,
rightmost_adj=rightmost_adj_set,
leading_node_index=leading_node
)
new_matrices = connect_adj_set(
leading_node_index=leading_node,
current_adj_mat=last_matrix,
adj_set=left_adj_set
)
for matrix in new_matrices:
if sum([sum(row) for row in matrix]) == sum_DS: # Check if the matrix is complete
complete_adj_mat.append(matrix)
else:
incomplete_adj_mat.append(matrix)
return complete_adj_mat
def show_graph_with_labels(adjacency_matrix, input_label=""):
"""
Displays a graph (network) represented by an adjacency matrix with optional labels.
Parameters:
- adjacency_matrix : the network structure.
- input_label : labels for the nodes of the graph (network).
"""
gr = nx.from_numpy_matrix(np.array(adjacency_matrix))
graph_pos = nx.spring_layout(gr, k=0.50, iterations=50)
# Draw nodes, edges, and labels based on the provided parameters.
nx.draw_networkx_nodes(gr, graph_pos, node_color='#1f78b4', node_size=220, alpha=0.6)
nx.draw_networkx_edges(gr, graph_pos, width=2, alpha=0.3)
if input_label:
labels = {i: str(label) for i, label in enumerate(input_label)}
nx.draw_networkx_labels(gr, graph_pos, labels)
else:
nx.draw_networkx_labels(gr, graph_pos)
plt.show()
def Network_score_with_cor(Candidate_network, Structure_score, Confidence_mat, Gene_list, Network_ranking_size, Permute_times = 10000, Permute_size = 2, Update_count = 1000):
"""
Calculates network scores based on network structure and correlation, applying permutations to get optimal gene labels.
Parameters:
- Candidate_network : potential network structures.
- Structure_score : Deviation measures for the structures.
- Confidence_mat : confidence matrix of data.
- Gene_list : genes of interest.
- Network_ranking_size : top n% potential networks for further analysis.
- Permute_times : maximum permutation times for each structure.
- Permute_size : number of genes to switch in each permutation.
- Update_count : threshold for increasing permutation size.
Returns:
- A tuple containing various outputs related to the network scoring and permutation process.
"""
tmp_structure_score = Structure_score.copy()
tmp_structure_score.sort()
final_gene_label = []
final_net_score = []
final_structure_list = []
Graph_list = []
unique_network_count = 0
structure_count = 0
structure_count_list = []
conv_net_score = []
update_net_score = []
Permute_size_list = []
Permute_gene_label = []
# execution
while unique_network_count < Network_ranking_size and structure_count < len(tmp_structure_score):
# ignore duplicated networks
net_index = np.where(Structure_score == tmp_structure_score[structure_count])[0][0]
structure_count += 1
current_net = np.array(Candidate_network[net_index])
# check isomorphism
isomorphic_result = False
if unique_network_count != 0:
#backward checking
current_graph = nx.from_numpy_matrix(current_net)
for iso_test in range(unique_network_count)[::-1]:
isomorphic_result = nx.is_isomorphic(Graph_list[iso_test], current_graph)
if isomorphic_result == True:
break
# calculate network score if network is unique
if isomorphic_result == False:
structure_count_list.append(structure_count)
Graph_list.append(nx.from_numpy_matrix(current_net))
final_structure_list.append(current_net)
unique_network_count += 1
current_degree_seq = sum(np.array(current_net))
# labeling
# start with max degree node
max_degree_node_index = np.where(current_degree_seq == max(current_degree_seq))[0][0]
# print('max degree node is located in "{}" index'.format(max_degree_node_index))
cor_sum = sum(Confidence_mat)
# print('sum of correlation for each gene: {}'.format(cor_sum))
max_cor_sum_gene = Gene_list[np.where(cor_sum == max(cor_sum))[0][0]]
# print('Labeling max degree node by gene ({}) with max sum of correlation.'.format(max_cor_sum_gene))
# (1) assign gene with max sum of correlation to node with max degree
gene_label = np.array(Gene_list.copy())
permute_candidate = np.array(np.where(gene_label == max_cor_sum_gene)[0][0])
permute_candidate = np.append(permute_candidate, max_degree_node_index)
# permute candidate one and candidate two
gene_label[permute_candidate] = gene_label[permute_candidate[::-1]]
# new gene co-expression matrix
gene_cor = np.copy(Confidence_mat)
gene_cor[:,permute_candidate] = gene_cor[:,permute_candidate[::-1]]
gene_cor[permute_candidate,:] = gene_cor[permute_candidate[::-1],:]
# (2) calculate initial score
net_score = 0
for jj in range(len(Confidence_mat)):
net_score = net_score + sum(gene_cor[jj]*current_net[jj])
# (3) permutation
# select two index to execute permutation
tmp_conv_net_score = []
tmp_Permute_size_list = []
tmp_update_net_score = []
tmp_Permute_gene_label = []
# Counter initialization to 0
no_update_count = 0
# Initialize total iterations
total_iterations = 0
tmp_permute_size = Permute_size
while total_iterations < Permute_times and tmp_permute_size < len(Gene_list):
tmp_gene_label = np.copy(gene_label)
permute_candidate = random.sample(range(len(Gene_list)), tmp_permute_size)
tmp_Permute_size_list.append(tmp_permute_size)
after_permute = permute_candidate.copy()
while after_permute == permute_candidate:
after_permute = random.sample(permute_candidate, tmp_permute_size)
# permute label
tmp_gene_label[permute_candidate] = tmp_gene_label[after_permute]
# new gene co-expression matrix
tmp_gene_cor = np.copy(gene_cor)
tmp_gene_cor[:, permute_candidate] = tmp_gene_cor[:, after_permute]
tmp_gene_cor[permute_candidate,:] = tmp_gene_cor[after_permute,:]
# calculate new score
tmp_score = 0
for jj in range(len(Confidence_mat)):
tmp_score = tmp_score + sum(tmp_gene_cor[jj]*current_net[jj])
# if score is higher than updating
if (tmp_score > net_score):
gene_label = tmp_gene_label.copy()
net_score = tmp_score.copy()
no_update_count = 0 # Reset counter to 0 after each update
tmp_Permute_gene_label.append(gene_label) # gene label after updating
else:
no_update_count += 1 # If no update, increment counter
tmp_conv_net_score.append(tmp_score)
tmp_update_net_score.append(net_score)
if no_update_count >= Update_count:
tmp_permute_size += 1
no_update_count = 0 # Reset counter to 0 when restarting the loop
i = 0
total_iterations += 1 # Increment total iterations
# (4) save result to a list
conv_net_score.append(tmp_conv_net_score)
update_net_score.append(tmp_update_net_score)
Permute_gene_label.append(tmp_Permute_gene_label)
final_gene_label.append(gene_label)
final_net_score.append(net_score)
Permute_size_list.append(tmp_Permute_size_list)
return final_structure_list, final_gene_label, final_net_score, conv_net_score, update_net_score, Permute_size_list, structure_count_list, Permute_gene_label
def truncated_power_law(alpha, maximum_value):
"""
Generates a discrete truncated power law distribution.
Parameters:
- alpha : the parameter of a power-law distribution.
- maximum_value : the maximum value of the distribution.
Returns:
- A sample of the distribution.
"""
x = np.arange(1, maximum_value + 1, dtype='float')
pmf = 1 / x**alpha
pmf /= pmf.sum()
return stats.rv_discrete(values=(range(1, maximum_value + 1), pmf))
